        SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1  RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
     2  TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,
     3  NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
     4  DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
C
        INCLUDE 'ABA_PARAM.INC'
C
        CHARACTER*80 MATERL
        DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1  DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2  STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3  PROPS(NPROPS),COORDS(3),DROT(3,3),
     4  DFGRD0(3,3),DFGRD1(3,3)
        double precision, DIMENSION(NTENS)::deviatoricStress
        double precision misesEqualStress,hydroStaticPressure,W
C
        DIMENSION EELAS(6),EPLAS(6),FLOW(6)
        PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
        DATA NEWTON,TOLER/1000,1.D-6/
        INTEGER K1,K2,cn

C
C -----------------------------------------------------------
C     UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC PLASTICITY
C     J2 FLOW THEORY
C     CAN NOT BE USED FOR PLANE STRESS
C -----------------------------------------------------------
C     PROPS(1) - E
C     PROPS(2) - NU
C     PROPS(3) - SYIELD
C     CALLS AHARD FOR CURVE OF SYIELD VS. PEEQ
C -----------------------------------------------------------

C
C     ELASTIC PROPERTIES
C
        EMOD=PROPS(1)
        EMU=PROPS(2)
        EG2=EMOD/(ONE+EMU)
        EG=EG2/TWO
        EG3=THREE*EG
        ELAM=EMOD*EMU/((ONE+EMU)*(ONE-TWO*EMU))
C
C     ELASTIC STIFFNESS
C
        DO  K1=1,NTENS
        DO  K2=1,NTENS
        DDSDDE(K2,K1)=0.0
        enddo
        enddo
C
        DO K1=1,NDI
        DO K2=1,NDI
        DDSDDE(K2,K1)=ELAM
        enddo
        DDSDDE(K1,K1)=EG2+ELAM
        enddo
        DO K1=NDI+1,NTENS
        DDSDDE(K1,K1)=EG
        enddo
C
C    CALCULATE STRESS FROM ELASTIC STRAINS
C
        DO K1=1,NTENS
        DO K2=1,NTENS
        STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
        enddo
        enddo
C
C    RECOVER ELASTIC AND PLASTIC STRAINS
C
        DO K1=1,NTENS
        EELAS(K1)=STATEV(K1)+DSTRAN(K1)
        EPLAS(K1)=STATEV(K1+NTENS)
        enddo
        EQPLAS=STATEV(1+2*NTENS)

        
C
C       HARDENING CURVE, GET YIELD STRESS
C
        CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS,NPROPS)

C
C       MISES STRESS
C
        call GetMisesEqualStress(STRESS,NTENS,NDI,misesEqualStress
     1  ,deviatoricStress,hydroStaticPressure)
C
C       DETERMINE IF ACTIVELY YIELDING
C
        IF (misesEqualStress.GT.(1.0+TOLER)*SYIEL0) THEN
C
C         FLOW DIRECTION
C       
        ONESY=ONE/misesEqualStress
        DO K1=1,NDI
            FLOW(K1)=ONESY*(STRESS(K1)-hydroStaticPressure)
        enddo
        DO K1=NDI+1,NTENS
            FLOW(K1)=STRESS(K1)*ONESY
        enddo
C
C       SOLVE FOR EQUIV STRESS, NEWTON ITERATION
C
        SYIELD=SYIEL0
        DEQPL=0.0
        RHS=misesEqualStress-EG3*DEQPL-SYIELD
        
        cn=0
        do while (ABS(RHS).GT.TOLER*SYIEL0)
            cn=cn+1
            DEQPL=DEQPL+RHS/(EG3+HARD)
            CALL AHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS,NPROPS)
            RHS=misesEqualStress-EG3*DEQPL-SYIELD
            if (cn.gt.10000) goto 10
        enddo
10      CONTINUE
C
C       CALC STRESS AND UPDATE STRAINS
C
        DO K1=1,NDI
            STRESS(K1)=FLOW(K1)*SYIELD+hydroStaticPressure
            EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO
            EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO
        enddo
        DO K1=NDI+1,NTENS
            STRESS(K1)=FLOW(K1)*SYIELD
            EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL
            EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL
        enddo
        EQPLAS=EQPLAS+DEQPL
        SPD=DEQPL*(SYIEL0+SYIELD)/TWO
C
C       JACOBIAN
C
        call GetMisesEqualStress(STRESS,NTENS,NDI,misesEqualStress
     1  ,deviatoricStress,hydroStaticPressure)
        W=9*EG**2/(misesEqualStress**2*(HARD+EG3))
        do i=1,NTENS
            do j=1,NTENS
                DDSDDE(i,j)=DDSDDE(i,j)-W*deviatoricStress(i)*deviatoricStress(j)
            enddo
        enddo     
        ENDIF
C
C STORE STRAINS IN STATE VARIABLE ARRAY
C
        DO K1=1,NTENS
        STATEV(K1)=EELAS(K1)
        STATEV(K1+NTENS)=EPLAS(K1)
        enddo
        STATEV(1+2*NTENS)=EQPLAS
C
        RETURN
        END
C
C
        SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,PROPS,NPROPS)
C
        INCLUDE 'ABA_PARAM.INC'
        DIMENSION PROPS(NPROPS)
        double precision A,B,C,E,SYIEL0
        PARAMETER (ONE=1.0D0)
C
C    SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
        E=PROPS(1)        
        A=PROPS(4)
        B=PROPS(5)
        C=PROPS(6)
        SYIEL0=PROPS(3)
C
C   IF MORE THAN ONE ENTRY, SEARCH TABLE
C
        if (EQPLAS.LE.0.0016) THEN
            HARD=1.0D8
            SYIELD=C+A*EQPLAS**B
        else
            HARD=A*B*EQPLAS**(B-ONE)
            SYIELD=C+A*EQPLAS**B
        endif
        
        RETURN
        END

        SUBROUTINE GetMisesEqualStress(STRESS,NTENS,NDI,misesEqualStress
     1   ,deviatoricStress,hydroStaticPressure)
        INCLUDE 'ABA_PARAM.INC'
        DIMENSION STRESS(NTENS)
        double precision, DIMENSION(NTENS)::deviatoricStress
        double precision misesEqualStress,hydroStaticPressure
        INTEGER NTENS,NDI

C     Calculate misese equal stress
        hydroStaticPressure=0
        do i=1,NDI
            hydroStaticPressure=hydroStaticPressure+STRESS(i)
        enddo 
        hydroStaticPressure=hydroStaticPressure/3
C     Calculate the deviatoric stress 
        do i=1,NDI
            deviatoricStress(i)=STRESS(i)-hydroStaticPressure
        enddo 
        do i=1+NDI,NTENS
            deviatoricStress(i)=STRESS(i)
        enddo 
        misesEqualStress=0
        do i=1,NDI
            misesEqualStress=misesEqualStress+deviatoricStress(i)**2
        enddo    
        do i=1+NDI,NTENS
            misesEqualStress=misesEqualStress+2*deviatoricStress(i)**2
        enddo  
        misesEqualStress=sqrt(1.5*misesEqualStress) 

        END      